Correlation Functions in SU(2)-Invariant Resonating- Valence-Bonds Spin Liquids on 

Nonbipartite Lattices 
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We introduce a Monte Carlo scheme based on sampling of PfafBans to investigate Anderson's 
resonating-valence-bond (RVB) spin liquid wave function on the kagome and the triangular lattice. 
This eliminates a sign problem that prevents utilization of the valence bond basis in Monte Carlo 
studies for non-bipartite lattices. Studying lattice sizes of up to 600 sites, we calculate singlet- 
singlet and spin-spin correlations, and demonstrate how the lattice symmetry is restored within each 
topological sector as the system size is increased. Our findings are consistent with the expectation 
that the nearest neighbor RVB states describe a topological spin liquid on these non-bipartite 
lattices. 



Introduction. It has been almost four decades since 
Anderson proposed[l| the quantum spin liquid state. Its 
undiminished appeal stems from a variety of applications 
from high temperature superconductivity 2] to quantum 
computing [3|, |j] . The nature of the short ranged variant 
of Anderson's "resonating valence bond" (RVB) spin liq- 
uid as a topological phase became understood through a 
series of papers [5|-M|. In particular, the invention of quan- 
tum dimer models [7| as an approximation to spin models 
finally lead to a lattice model exhibiting a topological 
RVB liquid phasejS], |24|. This however, did not immedi- 
ately address the (original) question whether this exotic 
phase could be stabilized within the phase diagram of 
S'?7(2)-invariant local spin-1/2 Hamiltonians. This was 
subsequently established for highly decorated lattices [9] 
and certain bipartite lattices [10|, by finding a parent 
Hamiltonian for the simplest, i.e., nearest neighbor ver- 
sion of the prototypical RVB spin liquid wave function on 
such lattices. Work on quantum dimer models|8l. llll - ll3l |. 
however, strongly suggests that nearest neighbor RVB 
states should be critical on bipartite lattices, as demon- 
strated recently [la, llll- They should describe a Z2-spin 
liquid with exponentially decaying correlations only in 
the nonbipartite case. While rigorously proven in the 
quantum dimer case, it is highly non-trivial establish this 
statement for the spin-1/2 RVB wave functions, due to 
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FIG. 1: a) Shape of the kagome lattice used in the calcula- 
tions. The lattice consists of m unit cells in the a direction 
and n unit cells in the b direction, for a total of 3mn sites. 
Periodic boundary conditions may or may not be introduced 
with periods ma and nb. b) The orientation used in the sign 
convention for the triangular lattice. 



local operators Oi and O2 takes on the form 
{RVB\0i02\RVB) Edd'{D\0i02\D') 



{RVB\RVB) 



J:d.d'{D\D') 



(1) 



orthogonality issues (cf, e.g., 
case, the nature of the correlation functions of the local 
spin and valence bond operator has not yet been studied 
systematically. This is largely due to a sign problem that 
will be addressed in this work. We finally mention that 
for the kagome case, the short-ranged RVB state studied 
here has been proven to be the ground state of a local par- 
ent Hamiltonian|18| (cf. also [19ll20|). The present work 
will provide essential evidence from correlations demon- 
strating that the kagome lattice RVB ground state of the 
Hamiltonian given in [18[ is a topological (Z2) spin liquid. 

Method. The standard method for calculating correla- 
tions of these wave functions on bipartite lattices is based 
on the observation that a general correlator between two 



Here D and D' represent dimerizations of the lattice, 
the sums run over all possible dimerizations, \D) is a 
nearest-neighbor valence bond (NNVB) state associated 
with a given dimerization and a link orientation of the 
lattice (defined below), and \RVB) is the RVB state, 
14|). In the nonbipartite \RVB) = J^d \^)- Since every pair of dimer config- 



D 

urations D,D' corresponds to a configuration of non- 
intersecting close packed loops on the lattice, Sutherland 
pointed out 2]J that the evaluation of such correlation 
functions may be reduced to the study of a classical loop 
gas model. This is so provided that the overlaps {D\D') 
are strictly non-negative. Otherwise, the evaluation of 
these correlators through Monte Carlo methods suffers 
from a sign problem. Indeed, it is not difficult to show 
that, e.g., for the kagome lattice, for any sign conven- 
tion for the states \D) some of the overlaps {D\D') are 
always negative. It is clear that such a sign prob- 

lem would never arise were we to work with an orthog- 
onal basis. In this case only strictly positive diagonal 
terms appear in the denominator of the expression re- 
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FIG. 2: The dimer-dimer correlation function shown for 
kagome lattices with PBCs and OBCs. Insensitivity to system 
size and short correlation length are evident. The PBC case 
has been calculated within a fixed topological sector. The 
inset shows a logarithmic plot including a linear fit, yielding 
a correlation length of 1.12(3). 



FIG. 4: The spin-spin correlation functions (SiSi+K.) and 
{SiS!+^) for difi'erent kagome lattices (PBC and OBC). 
Again, the topological sector was fixed in the PBC case. The 
inset shows a logarithmic plot with linear fit yielding a corre- 
lation length of 2.08(2). 



placing Eq. ([T]). An obvious candidate for such a ba- 
sis is the "Ising" -basis where local spins have definite z 
projection, \RVB) — J2i^i\^)^ ^^'^ ^ ''^^^ ^'^^r all pos- 
sible Ising spin configurations. Two primary questions 
need to be addressed to determine whether the Ising- 
representation lends itself to Monte Carlo evaluation of 
correlations. The first is the obvious question whether 
for the wave function \RVB), the coefficients a/ in the 
above representation can be efficiently calculated. The 
second question relates to the observation that for the 
short-ranged RVB state \RVB), it turns out that only a 
small fraction of configurations / will lead to non-zero aj. 
One may, however, ask if once an / with non-zero a/ is 
found, a sufficiently local update of / has a high chance 
of leading to a new /' with a/' ^ 0. To proceed, we first 
need to express the wave function \RVB) in the Ising ba- 
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FIG. 5: The expectation value of the dimer operator for 
links of the three possible directions and various lattice sizes. 
The average for one system size is shown as horizontal bar. A 
topological sector has been fixed. The discrepancy between 
inequivalent links rapidly decreases with system size, restoring 
the lattice symmetry. 
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FIG. 3: The dimer-dimer and spin-spin correlation func- 
tions for a 400 sites triangular lattice with OBCs. The inset 
shows logarithmic plot with fits, giving a correlation length of 
1.15(2) for the dimer-dimer decay. The spin-spin correlations 
display stronger even/odd effects at short distance. Fitting 
only odd distances in the spin-spin case gives a correlation 
length of 1.61(2). 



sis. We observe that in the Ising basis, the wave function 
\RVB) can be naturally written as a "Haffnian" [25J : 



ai = Haff[A/ij(/)] = 



1 



2^/2(f!) 



J2 M^,xM)Mx,xAl) 



xAh^_,xAl)- 



(2) 



AeSjs 



Here, M is a symmetric matrix whose indices run over 
the N lattice sites and which depends on the Ising con- 
figuration via Mij{I) = Qij{S„^^^S„.^i-Sa-i,i6„.^-f), where 
the ai describe the Ising configuration /. A runs over all 
permutations of the N sites, and Oy is the matrix de- 
scribing the chosen orientation of the lattice. An orienta- 
tion refers to a relation defined between any two nearest- 
neighbor lattice sites i, j, according to which either i < j 



or i > j holds, 
neighbors, 8.y = 



Then 8. 



if i, j are not nearest 



1 for i > j, and Oi 



-1 for i < j. 



Here we consider the orientation chosen for the kagome 
lattice indicated by the arrows in Fig. [T^. The formal 
definition of the Haffnian is related to that of the Pfaf- 
fian through omission of the sign factor (—1)'^. While the 
Pfaffian and the determinant can be evaluated in poly- 
nomial time, it is not known how to do this for the other 
two cases. It would thus be desirable to rewrite Eq. ([2]) 
through a Pfaffian. Luckily, this is the same problem that 
Kasteleyn solved long ago[22|, which has been a standard 
tool in the study of classical and quantum dimer mod- 
els. In the present context, it does not seem to have 
enjoyed much attention. Kasteleyn evaluated the parti- 
tion function of the classical dimer gas, which is exactly 
Eq. ([2]) with Mij replaced by |Gy|. He found that this 
problem may be written as PfafF[|8.y |8f^] where Q^ is a 
matrix similar to 0, but describing a different, so-called 
"Kasteleyn" orientation of the lattice. For planar lattice 
graphs, such an orientation may generally be found. A 
Kasteleyn orientation for the kagome lattice is given in 
[iat - The same method works in Eq. ([2]) [H. We thus 
have ai = Pia.S[Mij{I)Qfj]. We are now in a position to 
cast the problem of evaluating the correlation functions 
([1]) as a classical statistical mechanics problem. We have: 



{RVB\Oi0.j\RVB) 
{RVB\RVB) 



j:jj:j,ajaj,{I'\0,0,\I) 



(3) 

This may now be interpreted as the classical expectation 
value of a quantity /: (/) = ^^ fiG~^'/^i e~^' . Here, 
e~^' — |a/|^ and the value // of the quantity / in the 
Ising configuration / is given by // = J2i' i^'l^i^jl^)^- 

We have now demonstrated that the evaluation of cor- 
relation functions can be cast in terms of a partition func- 
tion, whose weights are positive and can be evaluated in 
polynomial time (the structure of our Pfaffian in fact 
allows reduction to the determinant of an N/2 x N/2 
matrix). Returning to our earlier caveat, we moreover 
found that once we have an initial Ising configuration 
/ with aj ^ 0, performing updates 25[ by exchanging 
neighboring spins has a high chance of leading to a new 
configuration /' with aj' ^ 0, The basic requirements for 
Monte Carlo evaluation are thus met. 

Results. Simulations are now performed for differnt 
lattices sizes. For the kagome lattice, we have chosen 
(m, n) as defined in Fig.[T]to be (10,5) for periodic bound- 
ary conditions (PBCs) and to be (20,8) and (20,10) for 
open boundary conditions (OBCs), resulting in a total 
number of iV = 150, 480, and 600 sites, respectively, 
(and in lattices with roughly unit perpendicular aspect 
ratio). For the triangular lattice, we show data belong- 
ing to a 20 X 20 "square" with diagonals (see Fig.[T]) giving 
a lattice of 400 sites. In one Monte Carlo sweep through 
the lattice, we attempt to do a number of N exchanges 



of two neighboring spins. All expectation values were 
calculated by making about 1,500,000 measurements on 
the configurations produced by the Monte Carlo process, 
allowing the system to equilibrate for 8000 sweeps. Au- 
tocorrelation times are generally quite low, on the order 
of 1. 

Fig. [5] presents the connected correlation function of 
the "dimer" or valence bond operator Si ■ 5^+^;, where 
i and i + x are nearest neighbors, for different lattice 
sizes and boundary conditions. It is evident that there 
is a finite and very short correlation length. From the 
inset it is clear that the absolute values of the correla- 
tion functions follow a simple exponential law already 
at short distance, from which we obtain a correlation 
length of ^ = 1.12(3). Moreover, the plot for 600 sites 
and OBCs coincides very well with that for 150 sites and 
PBCs. We note that for the case of PBCs, the method 
used to treat the classical dimer case{231 can again be 
adapted to the present situation, and yields an expres- 
sion of the amplitude a/ as a superposition of four Pfaffi- 
ans. Different such superpositions can be used to project 
onto different topological sectors of the toroidal system. 
While only one topological sector is shown, we have also 
convinced ourselves that results for different topological 
sectors agree within error bars. The fact that the dimcr- 
dimer correlations are apparently insensitive to both lat- 
tice size and boundary conditions, already for a relatively 
small size of 150 sites, is consistent with the hypothesis 
of a gapped state. We note moreover that the decay 
is very reminiscent of the quantum dimer model case, 
where dimer-dimer correlations have been shown to de- 
cay super-exponentially, with correlations being exactly 
zero beyond distance 2 26|. While this is clearly not the 
case for the RVB state, a very short correlation length 
of order 1 still mimics this behavior fairly closely. The 
qualitative agreement between the quantum dimer model 
and the RVB state is thus quite striking. 

Fig. [3] shows the dimer-dimer correlations for a 400 
site triangular lattice, displaying similarly short ranged 
correlations. Subdominant corrections to the dominant 
exponential decay are clearly somewhat more important 
than for the kagome, as one would generically expect; 
however a correlation length close to 1.6 is still clearly 
visible in the inset. All linear fits are obtained from a 
weighted least square regression, where the weights have 
been chosen as the inverse squares of the error bars. Note 
that although the value at distance zero has not been in- 
cluded into any fit, even this shortest distance data point 
tends to follow the exponential trend very well. We point 
out that a sign convention for the triangular lattice exists 
which eliminates the sign problem of Eq. ([T]) [27]. Here, 
however, we have chosen a different convention (Fig.[TlD), 
for which this problem persists. 

We also computed spin-spin correlation functions {Si ■ 
Sj). Results are shown for the kagome in Fig. 2] and 
for the triangular lattice in Fig. |31 Spin-spin corre- 



lations decay exponentially even in the critical square- 
lattice case[15|, and by theoretical prejudice should de- 
cay exponentially for all short ranged RVB states. More- 
over, even on the kagome, DMRG work has predicted 
a spin liquid phase with gapped spin but gapless sin- 
glet excitations 28[ . This might render the singlet sec- 
tor more crucial in the present context. Nonetheless, 
direct demonstration of the exponential decay of spin- 
spin corelations is not straightforward, especially in the 
presence of the sign problem discussed initially. Again, 
the short-ranged nature of the correlations is apparent in 
both cases. As a consistency check, both {Si ■ Sj) and 
3(S'fS'|) are shown, which must agree by SU{2) sym- 
metry. This symmetry is, however, not manifest in the 
Ising-basis we are working with. |36j 

Up to now we have demonstrated that connected corre- 
lations for the RVB states on the kagome and triangular 
lattice are short-ranged. This does, however, by itself not 
guarantee the liquid property of these states. In particu- 
lar, the four degenerate RVB ground states on the torus 
transform nontrivially under the space group of the lat- 
tice, and to demonstrate the liquid property and rule out 
the possibility of a valence bond solid [30| , it is essential 
to show that the full lattice symmetry is restored in the 
thermodynamic limit, for each individual ground state 
(within each topological sector) . We restrict ourselves to 
the kagome lattice here. In the following, we will refer to 
lattice links as "symmetry inequivalent" if they are not 
related by a symmetry of the wave function (even though 
they may be related by a symmetry of the lattice). For 
lattices of the shape shown in (Fig. [T^), with m, n both 
even, any three links along different directions will al- 
ways exhaust all possible classes of inequivalent links. In 
Fig. [5l we plot the expectation values of the dimer op- 
erator for three such links, evaluated in one topological 
sector, for various "even/even" lattices. One observes 
that the discrepancy between inequivalent links rapidly 
decreases, by a factor of at least 60 between 24 sites and 
48 sites, taking into account error bars. (The consistency 
between symmetry equivalent links suggests that the er- 
ror is much smaller than shown, and the factor is really 
on the order of 100). For larger lattice size, the calcula- 
tion becomes increasingly demanding, since, presumably, 
increasingly smaller error bars are needed to resolve the 
discrepancy in expectation values, while even maintain- 
ing the size of the error bars is more costly (Fig. [5|). It is 
worth noting, though, that the average of the three ex- 
pectation values for 72 and 96 sites appears to have con- 
verged, and we are thus approaching the thermodynamic 
limit. In all, these findings are highly consistent with 
the general expectation that the RVB-states describe a 
topological spin liquid. 

Conclusion. In this work, we have studied correlation 
functions of nearest neighbor resonating valence bond 
wave functions on both the kagome and the triangular 
lattice, with up to 600 lattice sites. A sign problem of 



earlier methods has been circumvented by using a Pfaf- 
fian representation of the wave function in the Ising basis. 
This allows for evaluation of correlators for both OBCs 
and PBCs, and, in the latter case, restriction to a single 
topological sector. This allowed us to present strong ev- 
idence that not only correlations decay exponentially as 
expected, but also that no broken lattice symmetry re- 
mains in the thermodynamic limit for the kagome lattice. 
For the kagome, this greatly adds to the amassed evi- 
dence that local SU{2) invariant Hamiltonians stabilizing 
a topological spin liquid state are possible[18|, |29[. Fur- 
ther possible applications of our method include the in- 
vestigation of short-ranged RVB wave functions on other 
nonbipartite lattices. In particular, certain next near- 
est neighbor links may be introduced in standard lattice 



geometries such as the square lattice[3l[, as long as the 



planarity of the lattice is maintained. This makes it nat- 
ural to introduce different weights for different types of 
valence bonds. Furthermore, our method allows for the 
introduction of any number of mobile (delocalized) holes 
and thus the study of monomer correlations and the re- 
lated confinement/deconfinement issue. We are hopeful 
that these prospects will stimulate future work. 

We would like to thank M. Ogilvie for insightful dis- 
cussions. This work has been supported by the Na- 
tional Science Foundation under NSF Grant No. DMR- 
0907793. Our MC codes are partially based upon the 
ALPS libraries [32il33|. 

Note added. After the completion of this work, Ref. [34] 
appeared using an approach that has been outlined in 
Rcf.[35]. 
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Supplemental material: Haffnian and Pfaffian 
representation of RVB state, and ergodicity 

We begin with the Haffnian representation of tfie RVB 
state \RVB) = Y.i<^Al) ■ We reproduce Eq. (2) of the 
paper as 



a/ = Haff[M,j(J)] 



1 



J2 Mx,xM)Mx,xAl) X • • • X AfA„_,A„(/) . (4) 

Independent of /, |A/y(/)| is just tlie adjacency matrix 
of the lattice, the only contributions come from permuta- 
tions A such that (A2n-iA2n) is a nearest neighbor pair, 
n being an integer running from 1 to half the number 
of lattice sites. Therefore, a contributing permutation 
A corresponds to a dimerization D of the lattice into 
nearest neighbor pairs. The permutations correspond- 
ing to the same dimerization in this way are all related 
by permuting the members of the individual pairs and by 



permuting pairs among each other. This over-counting is 
compensated by the combinatorial prefactor. The sum in 
Eq. (j4]) is thus essentially a sum over dimer configurations 
D. Since as function of cr\^^_^ and (Ta2„: -^A2„-iA2„(-^) 
is just the wave function of a singlet occupying the link 
(A2ri-iA2ri), cudowcd with the correct sign, the contribu- 
tion of a given A to Eq. ^ is just the spin wave function 
associated with the dimer covering D in the RVB state. 
On the other hand, the PfafRan representation can be 
written as 



a/ = Pfaff[My(/)] 



2^/2(f!) 



E i-^)'^lxAmlxAi) X ••• X ei^;_^,„(/) 



xes„ 



Mx,xM)Mx,xJI) X • • • X Ma„_,a„(/) . (5) 



For reasons of self-containedness, we reproduce here 
the well-known arguments[l| that Eqs. ^ and (O are 
interchangeable. This also serves to make it manifest 
that these arguments carry over from the original classi- 
cal dimer problem to the present situation without any 
difficulty. We first want to see that Eqs. (JU, ([5]) are 
identical within each topological sector of dimer cover- 
ings, possibly up to an overall negative sign. To this end, 
we use the fact that all coverings within a topological 
sector are mutually related by a sequence of "resonance 
moves" , where during a single resonance move, a set of 
dimers is shifted along a path of links forming a con- 
tractible (topologically trivial) loop. This interchanges 
occupied and unoccupied links along the loop. It is there- 
fore sufficient to show that the additional term of Eq. ([S]) 
compared to Eq. ^ is constant under a resonance move. 
Without any loss of generality, we can consider permu- 
tations A, A' corresponding to the dimer covering before 
and after the resonance move, respectively, that differ 
by a single cyclic permutation of the values of Ai . . . A2fc. 
Then (—1)^ = —(—1)'^ . Moreover, the defining property 
of the Kasteleyn orientation is that the number of clock- 
wise oriented links along any contractible loop of even 
length, which encloses an even number of lattice sites, 
is odd. (To ensure this, it is sufficient that the number 
of clockwise oriented links is odd around any elemen- 
tary plaquette of the planar lattice graph, for both even 
and odd plaquette perimeter). This then implies that 
the same cyclic permutation of indices described above 
also introduces an additional minus sign into the prod- 
uct Qxi.\2 ■ ■■■■ Qf2fc-i,A2fc- This shows that (g]) and ^ 
are identical within a topological sector up to a sign. 
For open geometry, there is just one topological sector 
(up to "non-resonable" configurations in the triangular 
lattice case [2] that are of negligible weight in the RVB 
state), and we are done. For periodic boundary condi- 
tions (torus geometry), Eqs. (JH) and (O do in general 
differ by a sign that is different for different topologi- 



cal sectors. However, this relative sign can be changed 
by further modifying Q along topologically non-trivial 
loops Q, making it possible to render all signs identical, 
or, as we did, to project onto a single topological sector 
by taking superpositions for different such modified Q^ . 
These considerations also set us up to observe an 
ergodicity property that the nearest neighbor spin ex- 
change moves employed in our Monte Carlo method have. 
Stated precisely, these moves are ergodic within the sub- 
set of Ising basis states that have non-zero weight for at 
least one nearest neighbor valence bond state (NNVB) 
within a given topological sector. First of all, it is ob- 
vious that nearest neighbor exchanges connect all Ising 
configurations belonging to a given NNVB state: Such 
configurations are linked by exchanges on links belong- 
ing to singlet pairs in the NNVB state. Moreover, for 
every dimer loop in the dimer covering associated with 
the NNVB state, there are Ising configurations contribut- 
ing to this state that are Neel ordered along this loop. 



Such Ising configurations, however, have non-zero weight 
not only in this particular NNVB state, but also in the 
NNVB state whose dimer configuration differs from the 
original one by a single resonance move along this partic- 
ular dimer loop. Then, by induction over the minimum 
number of resonance moves needed to connect two dimer 
coverings in the same topological sector, all Ising basis 
states having non-zero weight in the NNVB states asso- 
ciated with these coverings can be connected by nearest 
neighbor spin exchanges. 
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